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STEFAN'S  PROBLEM  IN  A  FINITE  DOMAIN  WITH  CONSTANT 
BOUNDARY  AND  INITIAL  CONDITIONS:  ANALYSIS 

Shunsuke  Takagi 
INTRODUCTION 

Under  any  initial  and  boundary  conditions  that  are  given  by  appro¬ 
priate  infinite  series,  Stefan's  problem  in  a  semi-infinite  domain  can  be 
solved  by  utilizing  the  two  infinite  series  of  the  elemental  temperature 
functions,  i.e. ,  a  general  solution  for  heat  conduction  in  a  semi-infinite 
domain  [1,  2,  3,  4],  Their  solution  method  was  brought  to  completion  in 
1978  by  Tao  [4],  He  introduced  a  formula  expressing  the  higher  derivatives 
of  a  function.  The  conditions  at  the  moving  boundary  can  be  rewritten  by 

this  formula  to  a  set  of  simultaneous  linear  equations  of  the  unknown 

parameters.  Further  progress  in  the  analysis  in  a  semi-infinite  domain  has 
been  made  by  Tao;  I  mention  a  few  that  may  be  interesting  to  us.  The 

density  jump  that  may  occur  at  the  moving  boundary  can  be  introduced  into 

the  analysis  (5).  Analyticity  of  the  interfacial  coordinate  as  a  function 
of  t^/^  proved  [6], 

In  spite  of  the  progress  that  has  been  made  so  far,  difficulties  in  a 
finite  domain  are  yet  outstanding.  Analysis  in  a  finite  domain  must  be 
understood  for  practical  applications.  The  heart  of  the  difficulties  seems 
to  be  in  analytically  formulating  the  well-established  practice  in  numeric¬ 
al  studies  (for  example  [7]),  i.e.,  that  the  solution  of  a  Stefan's 
problem  in  a  finite  domain  is  initially  close  to  that  in  a  semi-infinite 
domain  but  eventually  arrives  at  a  final  stationary  state.  In  this  paper, 
we  demonstrate  how  the  solution  in  a  finite  domain,  obtained  under  the 


simplest  boundary  and  initial  conditions,  transits  from  the  semi-infinite 
domain  solution  to  the  finite  domain  solution. 

The  breakthrough  is  achieved  by  use  of  three  essential  tools.  The 
first  is  the  conversion  of  a  pair  of  well-used  series  solutions  [1,  2,  3, 

4]  of  heat  conduction  in  a  semi-infinite  domain  to  a  pair  of  integral  ex¬ 
pressions.  One  of  the  expressions  enable  us  to  rewrite  Duhamel's  time  in¬ 
tegral  [8]  to  a  space  integral.  The  second  is  Widder's  [9]  integral  solu¬ 
tion  for  heat  conduction  in  a  finite  domain.  Its  use  enables  us  to  impose 
an  embedding  boundary  temperature  at  the  stationary  boundary  of  the  new 
phase  in  order  to  formulate  the  temperature  in  the  old  phase.  We  determine 
the  embedding  temperature  so  as  to  satisfy  the  interfacial  conditions.  The 
third  is  an  inverse-Laplace-integral  type  expression  of  inerfc  x  that  is 
valid  for  any  integer  n,  positive,  zero,  or  negative.  This  formula  is 
applicable  in  the  neighborhood  of  an  exponential  singularity. 

We  take  a  partial  sum  of  the  infinite  series  solution  for  the  tempera¬ 
ture  in  the  old  phase,  and  reinterpret  the  second  summand  contained  in  the 
last  term  of  the  partial  sum  in  the  following  way:  This  second  summand  is 
numerically  null  initially,  but  becomes  numerically  significant  at  an 
appropriate  time  and  continues  to  be  so  indefinitely  after  this  time.  If 
the  last  term  in  the  partial  sum  is  the  Nth  term  of  the  infinite  series,  we 
call  the  time  introduced  above  the  Nth  lead  time.  Prior  to  the  Oth  lead 
time,  the  first  term  approximates  the  solution  in  a  semi-infinite  domain. 

At  the  final  stage,  where  the  infinitely  many  lead  times  have  entered,  the 
temperatures  are  stationary  and  linear  in  both  the  new  and  old  phases, 
terminating  the  phase  change. 

The  first  part  of  this  paper  consists  of  the  mathematical  prelimi¬ 
naries  in  three  sections.  In  Section  1,  four  features  are  presented  with 
regard  to  the  elemental  functions.  First,  the  ranges  of  the  indexes  of  the 


elemental  functions  are  extended  to  negative  integers  so  that  the  formula 
for  the  series  development  of  a  function  of  a  series  may  be  smoothly 
applied  in  Section  3.  Second,  the  elemental  temperature  functions  are 
expressed  with  integrals  so  that  the  solution  for  heat  conduction  in  a 
semi-infinite  domain  may  be  given  an  integral  expression  in  Section'  2. 
Third,  an  inverse-Laplace-integral  type  expression  of  inerfc  x,  valid  for 
any  integer  n,  negative,  zero,  or  positive,  is  derived.  This  is  used  in 
Sections  7  and  9  for  the  transformation  in  the  neighborhood  of  the 
exponential  singularities  located  at  the  finite  terminal.  Fourth,  a  con¬ 
nection  with  Dirac's  delta  function  [10],  although  well-known,  is  presented 
so  that  it  may  be  smoothly  applied  in  this  paper. 

In  Section  2,  two  series  constituting  the  general  solution  for  heat 
conduction  in  a  semi-infinite  domain  [1,  2,  3,  4]  are  transformed  into  two 
integrals.  One  of  them  is  the  space-integral  expression  of  Duhamel's  time 
integral  [8,  p.  30]  for  solving  the  boundary-value  problem.  The  other  is 
the  well-known  integral  expressing  evolution  from  the  initial  value.  In 
Section  3,  a  formula  for  obtaining  a  serial  development  of  a  function  of  a 
series  is  presented.  The  elemental  functions  defined  on  the  moving  bound¬ 
ary  are  developed  into  a  series  of  t^^  by  this  formula.  The  formula  is 
applicable  at  a  nonsingular  point. 

The  second  part  of  this  paper  is  the  solution  of  Steran’s  problem  in  a 
finite  domain.  The  problem  we  solve  is  stated  in  Section  4.  In  Section  3, 
Widder's  [9]  solution  for  heat  conduction  in  a  finite  domain  is  transformed 
into  a  form  suitable  for  solving  the  problem  in  this  paper.  In  Section  6, 
the  embedding  boundary  temperature  with  unknown  coefficients  is  introduced 
at  the  stationary  boundary  of  the  new  phase  in  order  to  formulate  the 
temperature  in  the  encroached  old  phase.  In  Section  7,  the  infinite  series 
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formulating  the  unknown  temperature  in  the  old  phase  is  interpreted  as  an 
infinite  sequence.  Lead  times  are  introduced. 

In  Section  8,  the  boundary-value  problem  during  the  Nth  lead  time  is 
solved,  where  N  ^  0.  It  is  proved  that,  similarly  to  the  Neuman’s  solution 
[8,  p.  30],  each  of  the  solution  temperatures  consists  of  a  single  term. 
Prior  to  the  entrance  of  the  0th  lead  time,  this  solution  is  numerically 
equal  to  the  semi-infinite  domain  solution.  In  Section  9,  the  final  stage 
is  analyzed.  It  is  shown  that  the  temperatures  in  both  phases  are  station¬ 
ary  and  linear  in  space  at  the  final  stage.  Phase  change  therefore  finally 
stops. 


fATFE^ATICAI.  PRELIKINAPIFS 


1.  Elemental  Functions 

We  define  inerfc  x  for  a  nonnegative  integer  n  by 


irerfc  x  =  J 
/it  x 


(A-x)n 

n! 


f 


(1.1  ) 


where  -  00  <  x  <  ®  .  For  an  integer  0  V  £  n,  this  definition  yields 

d^  irerfc  x/dx^  =  (-1)^  ir  ^erfc  x  .  (1.2) 

Index  n  may,  therefore,  he  extended  to  negative  integers  by  defining 


i  nerfc  x  =  (-1)°  dnerfc  x/dx*1  .  (1.3) 

Then,  (1.2)  holds  true  for  any  positive  integer  k.  Using  the  Ferrite 
polynomial,  Fj,(x),  defined  by 


d 


k 


(-1) 


-x 


V*> 


(1.4) 


we  transform  (1.3)  to 


i  erfc  x  =  —  e  H  .  (x)  . 

r  n-1 

/it 


Tao  [4]  defines  C-n(x) 

G  (x)  =  {irerfc(-x)  +  (-l)n  iPerfc  x}  . 
n  z 

Substituting  from  (1.1),  we  transform  (1.6)  to 

00  2 

G  (x)  =  -----  /  (x-X  )P  e  dX  , 

T  nVm  -» 

which  integrates  to  an  nth  degree  polynomial, 


fn/Z]  (2x)n-2k 


_  t  ,  1  1  v  (2x) 

n  .,n  .  _r,  k !  (n- 


2  k=0 


2k)!  * 


The  derivatives  of  the  polynomials  Gn(x)  are: 


dkG  (x)/dxk  =  G  (x) 
n  n-k 


for  k  <  n 


for  k  *  n 


for  k  >  n 


Index  n  may,  therefore,  be  extended  to  negative  integers  by  defining 


G  (x)  =  0 
_n 


We  rewrite 


for  n  >  1  . 


k(x,<t)  =  e 


— X^  /  (4  ICt) 


//4lTKt  , 


to  an  inverse-Laplace-integral  type  expression, 

c+ 2 

k(x, <t)  =277  ■—  /  exp(x-  ~  ~l)dx  > 

/4<t  c-i00  /4  <t 

where  c  is  a  complex  number  whose  real  part  is  finite.  The  right  side  of 
(1.12)  transforms  to  the  right  side  of  (1.11),  when  the  former  is 
integrated  with  regard  to  p,  defined  by 


W  ■  (X  -  x//xt  )/2  . 

With  two  phases  in  Stefan's  problem,  we  use  Kt  in  place  of  t,  where  <  is 
the  thermal  diffusivity  of  a  phase  and  t  is  time. 

We  derive  the  formula 


c+i«  _  .2 

;rfc  x  =  — -  j  X  n  exp  (7 - Xx)dX  , 


(1.13) 


valid  for  any  integer  n,  -  ®  <  n  <  For  a  nonnegative  integer  n,  we  keep 
the  real  part  of  c  positive,  so  that  the  integral  is  convergent.  To  prove 
(1.13)  for  such  an  integer  n,  we  begin  by  transforming  (1.1)  to 


CO 

i°erfc  fx//4ictl  =  — r - - -  /  yn  k(x+y,  Kt)dy  . 

'  n!  (/4^F)n  0 


(1.14) 


Substituting  k(x+y,  <t)  from  (1.12)  and  evaluating  the  integral  with  regard 


to  y,  we  obtain 


.  c+i°° 

n  r  x  1  t 

erfc  -  =  —  X 

.  —  ■  TT1  J 


:)dx  , 


(1.15) 


which  is  equivalent  to  (1.13).  To  prove  (1.13)  for  a  negative  integer  n, 
we  begin  by  noting  that  (1.3)  is  equivalent  to 

i  nerf c( x//4<t )  =  2(-l)n  1  (/4<t)n  *  d°  *  k(x,<t)/dxn  * 

Substituting  k(x,Kt)  from  (1.12),  we  find  (1.15)  for  a  negative  n.  In  a 
finite  domain  U  ^  x  we  revise  (1.15)  to  a  nondimensional  form, 


/■/4Kt’\n  .n  ..  x 

l-j—J  1  erfc 

/4<t 


=  h  f  "  s'n_1  «XP ( £2  "  c|)dc 


(1.16) 


Equation  (1.12)  may  transform  to  an  alternative  form  [9,  p.  36], 


GO 

k(x,Kt)  =  -r—  J  exp(ixX  -  KtX^)  dX  . 


Letting  t  =  0,  the  above  becomes 


k(x,0)  =  ~  /  eixX  dX  . 

2  71  J 

—  oo 

Therefore  we  find 

Kx,0)  =  6(x)  , 

where  6(x)  is  Pirac's  delta  function  [10,  p.  38 ] . 

2.  Integral  Formulas 

The  heat  conation, 

3T/9t  =  <  32T/3y2  , 
has  a  general  solution, 

OO 

T(x,<t)  =  l  (/A  <t  )r  [A  G  (~)  +  B  ir  erf  c— — 1  . 

n=0  /4<t  /4  xt 

This  was  used  to  solve  Stefan’s  prohler  in  a  seri-infinite  domain  [1 

M. 

toe  ray  use  heat  functions, 
u^x.Kt)  =  n!(/4<t)n  i  nerfc  (>7/4<t ) 

and 

v^(x,Kt)  =  n!(/4<t)r'  (x//A<t  ) 

The  latter  is  called  the  heat  polynomial  hv  Kidder  [9,  p.  8-9). 
(1.14)  c,nd  (1.7),  (2.3)  and  (2.4)  rav  he  given  integral  expressions 
00 

v  (x,<t)  -  2  /  \,r  k(y+v,  xt  )dv 

n  f 

l; 

and 

OO 

f  F 

v^(x,xt)  =  J  y  Kx-v,  xt  )dy  . 


(1.17) 

(2.1) 

(2.2) 

,  2,  3, 

(2.3) 

(2.4) 

using 

(2.5) 

(2.6) 
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Therefore,  by  defining,  in  the  infinite  domain  -  00  <  x  <  «  ,  functions 


1 


i 

IS 

b 

-  . ' 

i 


f-z 

r.' 


f  (x) 


"  1 
-  I  b 


n=0 


A  x 
n !  n 


and 


E<X)  ’  2  l  h‘n*" 


n=0 


for  0  <  x  <  «> 


=  0  for  -  »  <  x  <  0 

(2.2)  transforms  to  an  integral  expression 


(2.7) 


(2.8a) 

(2.8b) 


T(x,Kt)  =  /  f(y)  k(x-y,  tct)dy  +  /  g(y)  k(x+y,  ict)dy  .  (2.9) 


Use  of  the  delta  function  (1.17)  shows  that  (2.9)  determines  the 
initial  value  f(x)  in  the  semi-infinite  domain  0  <  x  <  <»  ,  because  g(-x)  is 
null  for  x  in  this  domain.  The  initial  value  f(x)  in  the  integral  form 
(2.9)  may  be  discontinuous  and  is  more  general  than  the  serial  form  (2.2). 
We  assume  f(x)  to  be  odd, 

f (~x)  =  -  f ( x)  .  (2.10) 

Then  (2.9)  transforms  to 

00  00 

T( x, xt )  =  /  f(y)[k(x-y,  <t)  -  k(x+y,<t)]dy  +  /  g(y)  k(x+y,  <t)dy  .  (2.11) 

U  0 

The  first  integral  becomes  zero  at  x  =  0  for  any  t  in  the  domain  0  <  t  <  00 
and  the  second  integral  becomes  zero  at  t  =  0  for  any  x  in  the  domain  0  <  x 
<  «.  Therefore,  the  first  and  second  integrals  describe  the  evolutions 
from  the  initial  and  boundary  conditions,  respectively,  in  a  semi-infinite 
domain. 

Particularly,  we  have  an  expression  of  the  boundary  value, 

00 

T(0 ,  <t)  =  /  g(y)  k(y,tct)dy  .  (2.12) 

0 
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If 


<*>  _  n 

T(0,<t)  =  £  B  (/4xt)  inerfcO  ,  (2.13) 

n=0  n 

then  the  solution  g(y)  of  the  integral  equation  (2.12)  is  (2.8a). 

The  solution  of  the  boundary-value  problem  by  use  of  the  integral 
equation  (2.12)  is  consistent  with  Duhamel's  theorem  [8,  p.  30].  To  show 
this,  we  refer  to  Widder's  [9,  p.127]  solution  for  heat  conduction  in  a 
semi-infinite  domain.  In  our  notation,  it  is 


T(x, ict)  =  /  [k(x-y,  <t)  -  k(x+y,Kt)]  T(y,0)dy  + 
0 


t 

+  /  h(x,ic(t-x))  T(0,icx)dx  ,  (2.14) 

0 


where 

h( x, <t)  =  -2  3k(x,<t)/3x  =  (x/(<t))  k(x,  ict)  .  (2.13) 

The  first  integral  in  (2.14)  is  the  one  in  (2.11).  By  substituting 
T(0,<t)  from  (2.12),  the  second  integral  in  (2.14),  a  product  of  Duhamel's 
theorem,  becomes  the  repeated  integrals, 


«  t 

{  g(y)dy  /  h(x,<(t-x))  k(y , <x)d t  . 
0  0 


The  convolution  integral  in  the  above  simplifies  to, 


t 

/  h(x,tc(t-x))  k(y,xx)dx  =  k(x+y,*t)  ,  (2.16) 

0 


because  both  sides  yield  the  same  Laplace  transform, 


(l//4s”)  exp(-(x+y)//s ) 


To  derive  this,  we  employ  <t  in  place  of  t  in  the  usual  Laplace  trans¬ 
forms.  Widder's  solution  is,  therefore,  the  same  as  ours. 


3.  Functions  Of  A  Series 
At  the  interface. 


s(t) 


oo 


t(n+l)/2 


(3.1) 


the  temperatures  of  the  new  and  old  phases  need  to  be  expressed  as  func¬ 
tions  of  t^/2,  so  that  the  interfacial  conditions  may  be  expressed  as 
functions  of  t^^.  We  use  a  revised  version  of  Faa  de  Bruno's  formula 
[11,  p.  33]:  Let  functions  z  =  z(y)  be  analytic  and  y  =  y(x)  be  an  infinite 
series. 


OO 

(3.2) 

Then  the  power  series  expressing  the  composite  function  z(y(x))  is  given  by 


y-  l 

n=0 


z[y(x)] 


z[y(aQ)] 


where 


S  (a.)  -  T  ———————————— 

nv  j  [  X i!  ...  Xn_^l 


\i-vH 
an-  v+-l  » 


(3.3) 


(3.4) 


for  n  ^  1  and  z  v  1 ,  where  the  summation  is  over  all  the  sets  of 
at  least  one  nonzero  and  all  nonnegative  integers  Aj,  X2,  ...,  An  that 
simultaneously  satisfy  the  two  equations, 

Ai  +  X2  +  X3  +  ...  +  An-v+1  =  v  (3.5a) 

and 

A2  +  2A3  +  ...  +  (n-v)A  _  *  n  -  v  ,  (3.5b) 


for  every  integer  v  in  the  range 


1  <  v  <  n  .  (3.5c) 

The  solutions  are  tabulated  in  [12,  p.  831],  The  argument  a^  of  S^Ca^) 

stands  for  n  -  v  +  1  arguments  ai  , — , 

Substituting  s(t)  from  (3.1)  for  x  in  G  (x//4xt)  and  inerf c( x/ A Kt) , 

n 

we  find  the  series  of  t, 


(3.6a) 


and 


i'Wc  (-§^_)  -  l  l‘n>(4i-)’k  •  O-W 

/4xt  k=0  /4k 

The  coefficients  in  (3.6a)  and  (3.6b),  which  we  call  G-derivatives  and 
I-derivatives,  respectively,  are  given  for  k  =  0  by 


^(s./Ak)  =■  Gn(s0/Ax) 

and 


(3.7a) 


I{0n)(s./U<)  =  inerfc(sQ/Aic)  , 


(3.7b: 


and  for  k  >  1  by 


Gkn> 

K  /4k 


0 


I  g  h=-) 

v=l  n  v  Ax 


s.  (-dr) 


kv 


(3.8a) 


/4  k 


and 


i)n)  (-^r)  -  l  <->“  i""  Vfc  •  s.  h*-) 


S  . 


/4k  v=l 


/4k 


kv 


(3.8b) 


/4k 


The  index  n  -  v  can  be  negative.  Due  to  (1.10),  therefore,  the  number  of 
components  of  a  G-derivative,  shown  on  the  right  of  (3.8a),  may  be  less 


than  k  or  possibly  zero 


.  The  argument  s^/AIc  in  the  G-derivatives  and  the 

I-derivatives  stands  for  k  arguments  s^/Ak  , . s^//4<  . 

Substituting  x  in  (2.2)  with  s(t)  from  (3.1),  we  find  a  series  of  t. 


T(s(t),  Kt)  =  l  t[s)  tP  , 


(3.9a) 


where 


T(s)  =  l  (/47)n[A  G(n)  ^ — )  +  B  I(n)(-i-)j  . 
P  n=0  n  P_n  AT  n  P_n/4T 


(3.9b) 


Differentiating  (2.2)  with  x,  substituting  x  with  s(t)  from  (3.1),  we  find 
another  series  of  t. 


,  -  ft  (.(t),  «)  -  l  l'DS,TP 

3x  p=0  p 


(3.10a) 


where 


,(Ds)  _ 


1  f.  «(n-l )  /Sj 


l  (AT)  [A 

n=0  n  p  n  /4x  n  p  n  /4k 


(3.10b) 


Multiplication  of  T  by  3T/ 3x  on  the  left  side  of  (3.10a)  makes  the  series 
of  t  on  the  right  side  of  (3.10a)  conform  with  the  series  of  t  on  the  right 
side  of  (3.9a). 

We  call  (3.9b)  and  (3.10b)  level-p  coefficients  of  the  series  (3.9a) 
and  (3.10a),  respectively.  The  highest  index  of  the  parameters  in  the 
level-p  coefficients  are  p,  which  we  call  level-p  parameters.  Three  level- 


p  parameters  are  contained  in  the  level-p  coefficients.  Two  of  them  are  A 
and  Bp.  The  third  is  s^,  hidden  in  1^^  (s  j/Ak)  in  (3.9b)  and  in 
Ip  ^(sj//4k)  in  (3.10b),  as  may  be  found  by  use  of  the  formula, 


S  .(a  )  =  a  , 
nl  j  n 


(3.11) 


which  is  a  particular  case  of  (3.4).  In  the  equations  for  the  determina¬ 
tion  of  them,  which  we  discuss  subsequently,  the  higher-than-level-0  param- 


eters  are  all  linear  terms.  Only  sq,  one  of  the  level-0  parameters,  is 
nonlinear. 


STEFAN’S  PROBLEM  IN  A  FINITE  DOMAIN 

4.  Problem 

We  consider  the  simplest  Stefan's  problem  in  a  finite  domain  0  _<  x  _< 
l.  We  assume  freezing  starts  at  x  =  0  on  t  =  0.  The  boundary  temperature 
T^  at  x  =  0  is  a  constant  that  is  lower  than  the  freezing  temperature 
Tp  and  the  initial  temperature  Tg  is  a  constant  that  is  higher  than 
Tp.  The  boundary  temperature  at  x  =  l  is  a  constant  Tg.  We  have 

Ta  <  Tp  <  Tg  * 

At  t  =  0,  a  new  phase  emerges  at  x  =  0,  whose  temperature  we  express 
by  Tj(x,  Kjt),  where  tcj  is  the  thermal  diffusivity  of  the  new  phase. 

The  domain  of  the  new  phase  is  0  <  x  <  s(t),  where  s(t)  is  given  by  (3.1). 
The  temperature  of  the  old  phase  is  given  by  Tjj(x, xjjt) ,  where  kjj 
is  the  thermal  diffusivity  of  the  old  phase.  The  domain  of  the  old  phase 
is  s(t)  _<  x  l.  The  quantities  of  the  new  and  old  phases  are  designated 
by  the  roman  numerals  I  and  II,  respectively,  used  as  a  sub-  or  superindex, 
We  extend  the  domain 

s(t)  <  x  <  £ 

of  the  old  phase  at  time  t  to 

0  <x  <  Z 

by  introducing  the  embedding  boundary  temperature  at  x  =  0.  It  shall  be 
determined  so  as  to  satisfy  the  two  interfacial  conditions  at  the  freezing 
front  x  =  s(t). 


rr  f 

3T 

*1  ar  ■ 


9T_ 


3x 


Lp  dT 


(4.2) 


where  Kj  and  Kq  are  the  thermal  conductivities  of  the  new  and  old 
phases,  respectively,  p  is  the  assumedly  equal  density  of  the  two  phases, 
and  L  is  the  latent  heat. 

The  embedding  technique  was  initiated  by  Boley  [13],  Using  elemental 
functions,  ours  is  much  simpler  than  his. 


5,  Widder's  Solution  of  Heat  Conduction 

With  unrestricted  initial  and  boundary  conditions,  Widder's  [9,  p. 
131]  solution  of  heat  conduction  in  a  finite  domain  is  applicable  to 
solving  the  problem  of  embedding  posed  in  the  previous  section.  In  our 
notation,  his  solution  is 


Tn^x,Knt)  [e(x-y,Knt)  ~  eCx+y.K^t)]  T^Cy.o)  dy  + 


+  /  4*  ( x ,  <II ( t - T  ) }  Trr(0,KrrT)dT  + 


n 


(5.1) 


t 

+  /  ♦(*-X,Kn(t-T))  TnU,KnT)dT  , 


where 


0  (x ,  X£jt )  =  £  V(x  +  2n£,ic^t)  (5.2a) 

n=-°° 


and 


4>(x,ic^t)  =  £  b (x  +  2nfc,K^t)  .  (5.2b) 

n  =-°° 


Function  h(x,Kjjt)  is  given  in  (2.15).  We  call  the  three  integrals  on 
the  right  side  of  (5.1)  the  parts  due  to  initial  distribution,  embedding 


boundary,  and  terminal  boundary. 


The  first  integral  in  (5.1)  describes  the  evolution  fror  the  initial 


distribution  Tjj(x,  0).  To  show  this,  we  use  function  f(x),  which  is 
odd,  as  shown  in  (2.10),  and  is  periodic  vith  period  2£.  We  find  a  rela¬ 
tion  , 

00  £ 

/  f(y)  k(x-y,<IIt)dy  =  /  f(y)  [6(x-y,  ic^t)  -  0(x+y ,  ic  t )  ]dy  . 

-00  0 

Ve  assure  that  Tj-j(x,0)  is  extended  to  the  infinite  dorain  in  the  sane 
ranner  as  f(x)  is.  Letting  t  =  0  on  the  left  side  of  the  above  equation 
and  using  the  delta  function  (1.17),  and  noting  that  hy  letting  t  =  0  on 
the  right  side  of  (5.1)  the  second  and  third  integrals  both  disappear,  ve 
find  that  Tjj(x,  0)  is  the  initial  value. 

The  second  and  third  integrals  on  the  right  side  of  (5.1)  describe  the 
evolutions  fror  the  boundary  conditions  at  x  =  0  and  £,  respectively.  V;e 
clarify  this  staterent  in  two  parts:  The  forrer  or  latter  becores 
TH(0*Knt)  or  Tn(A»lfIIt)  lettlr£  X  *  0  or  t,  respectively. 

The  forrer  or  latter  reduces  to  zero  by  letting  x  =  £  or  0,  respectively. 
Note  that  the  first  integral  becores  0  at  x  =  P  and  x  =  £  because  of  its 
periodic  distribution. 

We  start  the  proof  by  forrulating  the  boundary  conditions  at  x  =  0  and 

A. 

00 

T  (0,K  t)  =  /  gj(y)  Ky,ie  t)dy  (5.3a) 

0  ■u' 

and 

00 

T  (£,k  t)  =  /  g2(y)  k(y,ic_t)dy  ,  (5.3b) 

0 

by  applying  (2.12)  to  two  seri-infinite  dorains,  0  <  x  <  •  and  £  >  x  >  -  °°, 
respectively,  where  y  in  (5.3a)  and  (5.3b)  stands  for  x  and  £-x,  respec¬ 
tively,  and  functions  gj(y)  and  g2(y)  are  defined  to  be  zero  for  negative 
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values  of  y.  The  delta  function  (1.17)  iray,  therefore,  be  applied  to  find 
the  Initial  terperature. 

We  substitute  (5.3a)  or  (5.3h)  into  the  second  or  third  integral, 
respectively,  on  the  right  side  of  (5.1),  divide  the  range  of  summation  on 
the  right  side  of  (5.2b)  into  the  nonnegative  integers  and  the  negative 
integers,  and  use  in  the  latter  the  oddity  of  the  h-function,  shown  hy 
(2.15),  so  that  (2.16),  in  which  the  arguments  rust  be  positive,  ray  be 
applied.  Thus  we  find 
t 

/  <J>(x,Kn(t-T))  Tn(0,<nT)dT 
00  00 

=  I  Sl(y)  I  {k(2n£+x+y,  k  t)  -  k(2(n+l )  £-x+y,  ic  t)  }dy  (5.4a 

0  n=0  11  LL 

and 

t 

/  Kn(t-T))  Tn(A,KI1T)dT 

00  00 

=  /  82<y)  1  {k((2n+l)£-x+y,<_t)  -  k((2n+l)  i+x+y,  ic  t  )dy  .  (5.4b) 

0  n=0  J-L 

If  we  let  x  =  0  or  l  on  the  right  side  of  (5.4a)  or  (5.4b),  we  find 
that,  after  the  cancellation,  only  the  n  =  0  term  remains,  which  is  equal 
to  T]j(0,icjjt)  by  (5.3a)  or  T^(  4,  *n;t)  by  (5.3b),  respectively. 

The  first  part  of  the  statement  is  thus  proved.  The  second  part 
immediately  follows. 

6.  Temperature  with  Embedding  Unknowns 

We  now  apply  our  initial  and  boundary  conditions  to  Widder's  unre¬ 
stricted  solution.  To  find  the  part  due  to  the  terminal  boundary,  we  let 
g2(y)  *  2  Tj  in  (5.4b)  and  use  (1.14)  with  n  **  0.  We  find  thus 


TB  =  T_  l  [erf 


(2n+l ) £-x  ,  (2n+l)£+* 


/4<nt 


-  erfc^LZL^J  . 


(6.1a) 


A<nt 


To  integrate  the  part  due  to  the  initial  distribution,  we  let  Tjj(y,0)  = 

Tg  in  the  first  integral  on  the  right  side  of  (5.1),  substitute  the 
0-functions  from  (5.2a),  change  the  range  of  integration  from  0  to  £  to  the 
difference  of  the  one  from  0  to  »  and  the  one  from  £  to  “  ,  carry  out  the 
thus  defined  integrations  by  use  of  (1.14)  with  n  =  0,  and  change  those 
with  negative  argument  to  those  with  positive  argument  by  use  of  (1.6).  We 
find  thus 


TP  =  T  [1  -  I  (-l)Perfc  ^  (_i)%rfc  . 

n=0  /4<j_jt  n=l  /4  t 


(6.  Ih) 


Adding  (6.1a)  and  (6.1b),  we  get 


TB  +  ID  =  Tr{  1  -  l  erfc  +  £  erfc  ----- ' 

n=0  Mcjjt  n  =  l  /4<^t 


(6.1c) 


I'sirg  g(y)  in  (2.8)  for  g^fy)  in  (5.4a),  rewriting  n  and  Bn  in  tbe 

n  V 

initial  condition  (2.13)  to  k  and  P^  /£  ,  respectively,  and  using  (1.14), 

w’e  integrate  tbe  part  due  to  tbe  embedding  boundary  condition. 

Summing  tbe  results,  we  obtain 

*  T„  + 


(/ 4  K  t  V  00 

Jo 


k  r  2n£+x  .V  2(n+l)£-xi 

erfc  -  -  i  erfc - - — 


(6.2) 


/4<nt 


A<ut 


n  n 

Tbe  coefficients  ,  Pj  are  to  be  determined.  Included  in  tbe  sur- 

mand  for  k  =  0  in  (6.2),  tbe  forrula  in  tbe  pair  of  braces  in  (6.1c)  need 
not  be  listed.  It  is  easy  to  directly  check  that  (6.2)  satisfies  all  tbe 
conditions  thus  far  imposed.  Differentiation  yields 
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9T, 


3x 


(x,Knt)  - 


00  B.  //4  k  t  V  «°  ,  ,  _  . , 

v  V  I  U  1  r  r  . ^ —  1  2n£+y 

l  7=  \—£ — /  l  i1  erfo  n= 

k=C  /4k^  '  n=0  /4<^t 


,  . V-l  ,  2(p  +  ] ) £-x 1 

+  1  erfc — -  — } 

A*nt 


(6.3) 


7 .  I ead  Tires 

We  rev  rite  (6.2)  to  a  secuence 
Tn(x,  Knt)  -  lir  ^’(x,  K  t)  , 

N-voo 


(7.1) 


where 


Tn)(x’Kn° 


tb  + 


00  /4ic  t  k  N 

+  l  B™  (-T3-)  I  |ikerfc  -  lkerfc2(^‘  >*-*■)  , 


k=0 


n=0 


/4Knt 


/4xnt 


(7.2) 


which  satisfies  the  terminal  boundary  condition  and  the  initial  condition. 
UN 

If  coefficients  are  chosen  to  satisfy  the  interfacial  conditions 

(4.2),  we  call  (7.2)  the  Nth  stage  solution. 

To  understand  the  successive  emergence  of  infinitely  many  stages  in 

(7.1),  we  introduce  a  lead  time  in  the  Nth  stage  solution.  We  define  it  in 

UN 

a  special  case  such  that  B^  =0  for  all  the  integers  k  >_  K  +  1  .  Time  t 


that  makes  i  erfc [ (2(N+1 ) l  -  s(t)  )/ /4<^t  ]  numerically  effective  will  be 
called  the  Nth  lead  time.  Note  that  this  is  the  least  of  the  values  that 
i  erfc[  (2(N+1  )£-x)//4<^t  ]  in  the  last  pair  of  braces  in  (7.2)  can  take  for 
k  =  K  in  the  domain  s(t)  <  x  <  I  at  time  t. 

This  definition  is  appropriate  to  the  problem  we  are  solving  in  this 

UN 

paper,  where,  as  we  shall  show,  all  the  coefficients  B^  are  zero  for  k  _> 

UN 

1  and  therefore  only  Bq  is  nonzero.  A  lead  time  is  a  semi-infinite  time 


interval.  Prior  to  the  0th  lead  time  in  our  problem,  the  0th  stage  solu- 
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tion  is  the  semi-infinite  domain  solution.  When  the  Oth  lead  time  is 
effective,  the  second  summand  exists  in  the  0th  stage.  The  final  stage  is 
arrived  at  in  the  infinite  stage. 

Although  not  a  constant  but  a  time  interval,  a  lead  time  nust  be 
dealt  with  like  a  constant  in  the  time  differentiation.  Except  prior  to 
the  Oth  lead  time  in  the  Oth  stage,  the  transcendental  equation  for  the 
determination  of  sQ  at  the  Nth  stage  includes  the  lead  time,  as  we  shall 
show,  and  therefore  the  root  s0  is  a  function  of  a  lead  time.  We  express 
s(t)  at  the  Nth  stage  by 


N 

Y 

L 

n=0 


where  coefficients 
Similarly  to  (7.2) 


,  we 


c(n+1>/2  ,  (7.3) 

are  in  general  functions  of  the  Nth  lead  time, 
rewrite  (6.3)  to  the  expression  at  the  Nth  stage. 


"S’ 

T  —  <x-'!nt)  ■ 


ttn  - 

=  -Y  Y  {ik_1erfc  +  ik~1erfc^S-+1)— xl.  (7.4) 

k=0  AkZ-  £  n=0  SHZt  AZT? 


•n  "  -  Knt  fHV 

To  find  the  series  of  x  at  the  interface,  we  use  the  formula  (1.13)  to 
obtain  the  integral  expression 

2nl 


2n=  ±  n)  -  ^  rv-1  exp(£  -  x  _ 

/4<_t  c-i-  ^S-rt 


i  erfc(- 


±  Xn)dX  ,  (7.5) 


where 


c  1  (N)  p 

n  -  l  - s_  x 


P=0  p 


Use  of  the  development  (3.3)  yields  the  expansion 

(N)\  /_<N) 


+  Xn 
e  “  exr 


r+Xs 


0 


OD  n 

I  TP  l  (* 


x)v  s 


li 


(4<^  /  p=0  v»0 


pvl 


(7.6) 


/4t 
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where  the  convention 


WaJ>  ' 1 

and 

S  (a. )  =  0  for  n  >  1  , 

nO  j  — 

is  used.  Substituting  (7.6)  in  (7.5)  and  using  (1.13),  we  obtain  the 
series  development, 


.k  ^  2n£  ±  s^N\t)  v 

l  erf c  -  =  l 


/4<nt 


p=0 


l  (+D 


v=0 


„  2„E  i  ,<« 

.  k-v  0 

l  erf c - 


✓4 


Kn' 


(N) 

Letting  x  =  s  (t)  in  (7.2)  and  (7.4)  and  using  (7.7),  we  find 


of  x: 


Tn<s(N,(t),  «nt>  ■  I.®'’1  ,P 


and 


3T__  « 

x  •  (s(  (t),  K—t)  »  i  T 


3x 


p=0 


(HN,Ds)  p 
P 


We  neea  only  the  first  terms. 


TmN,s>  ,  Tb  t 

*  B°N  J  {erfc 
n=0 


2nt  ♦  s‘N)/7 


erfc- 


2(n+l)£  -  s^Vt 


/4<nt 


/4xnt 


and 

( ElN.Ds)  _ 

0 


HN 


N  2n£  +  sl^Vt  2 ( n+ 1 ) £.  - 

Y  {i-1erfc - - -  +  rlerfc - - -°  )  . 


/t 


(7.7) 

series 


(7.8a) 


(7.8b) 


(7.9a) 


(7.9b) 


Because  the  coefficients  for  p  _>  1  are,  as  we  shall  show,  all  equal  to 

zero,  the  higher  terms  need  not  be  made  explicit. 


8.  The  Nth  Stage 

The  general  solution  (2.2)  gives  the  Nth  stage  temperature  of  the  new 
phase. 


T^^x.Kjt)  =  l  (AlCjt)  [A™  Gr(-^)  +  irerfc(~=)]  .  (8.1) 

n  =0  /  4  k  _  t  /4  x  t 


IN  .  n 


Substituting  x  with  s  (t)  from  (7.3),  and  expanding  the  elemental  func- 

(N)  (N) 

tions  by  use  of  (3.6a)  and  (3.6b),  we  rewrite  T  (s  ( t ) ,  x^t)  to  a 


series  of  x, 


V)  -  l 


(8.2a) 


where 


I  A™  G(">/-1l\+  .  (8.2„, 

k  „io  -  k-"\//nrl  ”  >-r\Ar/ 


Differentiating  (8.1)  with  x,  substituting  x  with  s^  ;(t)  from  (7.3),  and 


using  (3.6a)  and  (3.6b),  we  find 


•^(s<N)«>,  V)  .  z 


(8.3a) 


where 


T<IN'Ds)  -  I  (/kX-)"'1  AIf'  C 
k  Ln  K  1 J  n 


l"n  \/xrJ  ”  l"n 


(8.3b) 


Rewriting  the  boundary  conditions, 

tJN)(0,  xTt)  =  T 


(8.4a) 


T{N)(s(N)(t),  <  t)  =  T 


(8.41 


N  N  ( 1 

into  series  of  t,  we  can  successively  determine  A  and  B  as  functions  s' 

p  r  o 

(N) 

s  .  The  first  pair  of  this  seouence  is: 

P 


=  rF  -  ta  erfc^soN)//4,ci^rf (sJK)//4ki) 


(8.5a 


B0  =  ~  (TF  "  * 


(8.5b 


Ve  write  Erf  vdth  capital  F  lest  it  he  confused  with  erfc  x.  The  higher 
pairs  of  the  sequence  need  not  be  described  in  our  prohler,  where  they  are 
as  we  shall  show,  all  equal  to  zero. 


Substituting  T 


(IIN.S) 


froir  (7.9a),  the  coefficient  of  t  in  equation 


(4.2a)  yields  the  evaluation  of  P 


■  -  <TB  -  V'V 


(8.6a 


where 


N  2n£  +  s„  /t 

Rj,  =  J  [erfc - —  — -  -  erfc 

n=0  /4<_t 


2(n  +  l )  i  -  sfVt 


/4Kn£ 


(8.6b 


Substituting  found  fror  (8.3h)  with  k  =  0  and 

from  (7.9b),  the  coefficient  of  in  equation  (4.2b)  yields  the  transcen- 

(N) 

dental  equation  for  the  determination  of  sn  , 


where 


qn  •  I  t1'1"* 

n=0 


2ni  +  sn  /t  2(n+l)£  -  s~  Vt 

: - - g: - +  i^erfc - 2 - ]  . 

/  4  <  ^  t  /4  t 


(8.8) 


(N)  IN 

Thus  determined,  Sq  is  a  function  of  the  Nth  lead  time.  Then,  Aq  , 

IN  UN 

Bq  ,  and  Bq  are  also  functions  of  the  Nth  lead  time. 

We  have  thus  evaluated  the  level-0  parameters.  We  call  the  equation 
for  determining  the  level-n  parameters  level-n  equations.  For  n  1,  the 
level-n  parameters  are  linear  in  the  level-n  equations.  The  level-L  equa¬ 
tions  are  homogeneous  with  respect  to  the  level-1  parameters,  which  are 
therefore  all  equal  to  zero.  Proceeding  this  way,  we  find  that  the  higher- 
than-level-0  parameters  are  all  equal  to  zero. 

Thus  we  find 


TjN)(x,Klt)  -  +  (I,  -  T  )  Erf  -jb/ 

/4<_t  ' 


(8.9a) 


I  ,  .x  1  „IN  r,_i  ,  x  ] 

T  •  ^  (x,Kit)  --  —  B0  [i  erfc  —  ] 

/  4<j  /4iCjt 


(8.9b) 


for  the  new  phase.  Also  we  find 


4N><x,*nt)  -  Ts  +  B™  X  [erfc  -  erfc^j^] 

n=0  /4  Kjjt  /4icnt 


(8.10a) 


3Tn 

3x  (x*’ent) 


for  the  old  phase. 


/  4  n=0 


l  [i^erfc  +  i-W2-^^) 


/4<nt 


/4Knt 


(8.10b) 


9.  The  Final  Stage 

To  formulate  the  final  stage,  we  transform  the  infinite  series, 


,  r/4<t  r  r.k  „  2n £+x  .k  2(n+l)£-xi 

Lk  =  l — £ — J  1  U  erf e- — 1 "  -  i  erfc— — J 


(9.1) 


w  <V4Kt  -*k  r  r,k  2n£+x  .  ,k  r  2(n+l)£-x  i 

\  =  (  £  )  I  [i  erfC"=^  +  1  erfe~P= - 1  * 

n=0  /4<t  /4xt 


(9.2) 


contained  in  (6.2)  and  (6.3),  respectively.  We  drop  the  subscript  H  in 
this  computation.  Although  our  problem  is  concerned  with  Lq  and  M_j,  we 
transform  the  general  cases  and  at  this  stage  of  the  exposition. 

We  begin  with  the  transformation  of  the  infinite  series. 


<V4ict-ik  r  .k  ..  2n£+x 

l — T~ )  l  i  erfc  -^7- 


Using  formula  (1.16),  this  becomes 


.  c+i®  .  .  « 

1  ,  -k-1  ,  ,  Kt  x^  _ 

-  7  / «  ^<{2  7  - « 7)  L 

c-i®  £  n=0 


C+i® 

1  ,  -k-1 


‘2nl  /  5  exp( C2  —  +  £  t  )  Cosechc  d£ 

c-i®  * 


We  write  Cosech  x  with  capital  C,  lest  it  be  confused  with  cosec  x. 
Similarly, 


,/4xt.k  „  k  2(n+l)£-x 

(  — )  I  i  erfc  — — - 

n“0  /4  <t 


transforms  to 


1  C+l~  -k-1 


/  C  exp  (Cz  —  -  ^  t  }  Cosech  £  d£ 


^  £2 


24 


Thus  we  find 


i  c'fl*  -k-i 


5 

c-i® 


Kt 


£-x 


exp(c2  — )  Sinh(?  — “)  Cosech?  dE,  . 
I2  % 


(9.3a) 


Expanding  the  Sinh  function  by  use  of  the  addition  formula,  and  using 
(1.16)  and  (1.6),  we  may  give  (9.3a)  another  form. 


L.  = 


\/4<t  / 


k  _ x_ 

i  erfc 


/4  (ct 


c+i® 


-  ^  J  fkl  exp(c2  Sinh(-f£)  CothC  dC 


(9.3b) 


c-i® 


V 


When  x  is  substituted  for  s(t)  from  (3.1),  the  latter  transforms  to  a 
series  of  t,  but  the  former  does  not.  Similarly,  we  find 


i  ,i+1*  -k-i 


Kt  . 


l-x. 


^  =  iri  /  E,  exp(c2  — )  Cosh(C  l  )  CosechE  d£ 


(9.4a) 


c-i® 


and 


V4xt 

\  =  -  (  — ) 


/  x  \  l+(-Q* 

c\ /Kt  /  "  2 


k  _x 
i  erfc- 


AZ 


l  -k-i 


Kt  , 


+  ~  J  C  “  *  exp(cz  ^-)  Cothc  Cosh  4s  <U  . 


iri 


(9.4b) 


c-i®  l 

In  our  problem,  the  temperatures  at  the  final  stage  are  stationary  and 
linear  in  space.  To  show  this,  we  begin  by  letting 


)(£  l 

=  exp( £Z  — )  Sinh(c  ~ ySinh  E, 


(9.5a) 


in  (9.3a)  and 
<KS)  *  exp 


/  0  Kt  \  /  i~X 1 

U2  — ^  )  Cosh  ( 5  t  1( 


C/SinhC) 


(9.5b) 


in  (9.4a).  Noting  that  the  $(5)’s  are  even  functions,  we  find  that  I<) 


and  M-i  are  given  by  the  respective  residues  at  £  =  0.  Thus  we  find 

Lq  =  1  -  x/ l  (9.6a) 


M_i  =  1  . 

To  formulate  the  final  temperature,  we  first  note  that 

lim  s<N)  =  0 

Nx» 

and  that  the  final  interface  coordinate,  s(“),  is  given  by 


s(<*>)  =  lim  S<,N)  /t 

N-h» 


(9.6b) 


(9.7) 


Note  that  as  N  -*■  <*>  y  also  t  <*>.  Taking  the  limit  in  (8.9a),  and  noting 


Erf(x)  '=  —  x  , 


(9.8) 


for  a  small  x,  we  find  that  the  final  temperature  in  the  new  phase  is  given 


Tj0°) (x, Kjt)  =  Ta  +  (Tf  -  Ta)x/s(.)  . 

Using  M_i  in  (9.6b)  and  Lq  in  (9.6a),  we  find 


(9.9) 


lim  Q„  - 
N-h» 


/4Knt/£ 


(9.10a) 


lim  Rjj  ■  1  —  s(®)/£. 

N+oo 


(9.10b) 


Dividing  all  the  terms  in  (8.7)  hy  t1'  ,  and  letting  t  ♦  °°,  ve  find  that 
the  limit  of  (8.7)  is  stationary,  given  by 

(Tb  -  T 


The  phase  change  therefore  finally  stops.  Using  (9.10b)  in  (8.6a),  we  find 
that  the  limit  of  (8.10a)  is  given  by 

T^Oc^t)  =  Tr  -  (Tb  -  -  s( °°) ) .  (9.12) 

CONCLUSIONS 

Py  using  the  simplest  boundary  and  initial  conditions  we  have  shown 
that  Stefan's  problem  in  a  finite  domain  can  be  solved.  To  solve  the 
probleir,  we  have  erployed  the  solution  for  heat  conduction  in  a  finite 
doirain  whose  elemental  temperature  functions  are  members  of  the  family  of 
the  error  function.  The  effect  of  the  finite  terminal  is  the  introduction 
of  singularities.  Starting  at  a  solution  in  the  semi-infinite  domain  and 
passing  through  infinitely  many  lead  times,  the  solution  in  the  finite 
domain  arrives  at  the  final  stationary  stage. 
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